/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef __itkDifferenceImageFilter_hxx
#define __itkDifferenceImageFilter_hxx

#include "LOCAL_itkDifferenceImageFilter.h"

#include "itkConstNeighborhoodIterator.h"
#include "itkImageRegionIterator.h"
#include "itkNeighborhoodAlgorithm.h"
#include "itkProgressReporter.h"

namespace itk
{
// ----------------------------------------------------------------------------
template <class TInputImage, class TOutputImage>
LOCAL_DifferenceImageFilter<TInputImage, TOutputImage>
::LOCAL_DifferenceImageFilter() :
  m_DifferenceThreshold(NumericTraits<OutputPixelType>::Zero),
  m_MeanDifference(NumericTraits<RealType>::Zero),
  m_TotalDifference(NumericTraits<AccumulateType>::Zero),
  m_NumberOfPixelsWithDifferences(0),
  m_ToleranceRadius(0),
  m_IgnoreBoundaryPixels(false)
{
  // We require two inputs to execute.
  this->SetNumberOfRequiredInputs(2);
}

// ----------------------------------------------------------------------------
template <class TInputImage, class TOutputImage>
void
LOCAL_DifferenceImageFilter<TInputImage, TOutputImage>
::PrintSelf(std::ostream & os, Indent indent) const
{
  this->Superclass::PrintSelf(os, indent);
  os << indent << "ToleranceRadius: " << m_ToleranceRadius << "\n";
  os << indent << "DifferenceThreshold: " << m_DifferenceThreshold << "\n";
  os << indent << "MeanDifference: " << m_MeanDifference << "\n";
  os << indent << "TotalDifference: " << m_TotalDifference << "\n";
  os << indent << "NumberOfPixelsWithDifferences: "
     << m_NumberOfPixelsWithDifferences << "\n";
  os << indent << "IgnoreBoundaryPixels: "
     << m_IgnoreBoundaryPixels << "\n";
}

// ----------------------------------------------------------------------------
template <class TInputImage, class TOutputImage>
void
LOCAL_DifferenceImageFilter<TInputImage, TOutputImage>
::SetValidInput(const InputImageType *validImage)
{
  // The valid image should be input 0.
  this->SetInput(0, validImage);
}

// ----------------------------------------------------------------------------
template <class TInputImage, class TOutputImage>
void
LOCAL_DifferenceImageFilter<TInputImage, TOutputImage>
::SetTestInput(const InputImageType *testImage)
{
  // The test image should be input 1.
  this->SetInput(1, testImage);
}

// ----------------------------------------------------------------------------
template <class TInputImage, class TOutputImage>
void
LOCAL_DifferenceImageFilter<TInputImage, TOutputImage>
::BeforeThreadedGenerateData()
{
  ThreadIdType numberOfThreads = this->GetNumberOfThreads();

  // Initialize statistics about difference image.
  m_MeanDifference = NumericTraits<RealType>::Zero;
  m_TotalDifference = NumericTraits<AccumulateType>::Zero;
  m_NumberOfPixelsWithDifferences = 0;

  // Resize the thread temporaries
  m_ThreadDifferenceSum.SetSize(numberOfThreads);
  m_ThreadNumberOfPixels.SetSize(numberOfThreads);

  // Initialize the temporaries
  m_ThreadDifferenceSum.Fill(NumericTraits<AccumulateType>::Zero);
  m_ThreadNumberOfPixels.Fill(0);
}

// ----------------------------------------------------------------------------
template <class TInputImage, class TOutputImage>
void
LOCAL_DifferenceImageFilter<TInputImage, TOutputImage>
::ThreadedGenerateData(const OutputImageRegionType & threadRegion, ThreadIdType threadId)
{
  typedef ConstNeighborhoodIterator<InputImageType>                           SmartIterator;
  typedef ImageRegionConstIterator<InputImageType>                            InputIterator;
  typedef ImageRegionIterator<OutputImageType>                                OutputIterator;
  typedef NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> FacesCalculator;
  typedef typename FacesCalculator::RadiusType                                RadiusType;
  typedef typename FacesCalculator::FaceListType                              FaceListType;
  typedef typename FaceListType::iterator                                     FaceListIterator;
  typedef typename InputImageType::PixelType                                  InputPixelType;

  // Prepare standard boundary condition.
  ZeroFluxNeumannBoundaryCondition<InputImageType> nbc;

  // Get a pointer to each image.
  const InputImageType *validImage = this->GetInput(0);
  const InputImageType *testImage = this->GetInput(1);
  OutputImageType *     outputPtr = this->GetOutput();

  if( validImage->GetBufferedRegion() != testImage->GetBufferedRegion() )
    {
    itkExceptionMacro( << "Input images have different Buffered Resions." )
    }

  // Create a radius of pixels.
  RadiusType         radius;
  const unsigned int minVoxelsNeeded = m_ToleranceRadius * 2 + 1;
  const typename TInputImage::SizeType imageSize = validImage->GetBufferedRegion().GetSize();
  for( unsigned int d = 0; d < TInputImage::ImageDimension; ++d )
    {
    if( minVoxelsNeeded < imageSize[d] )
      {
      radius[d] = m_ToleranceRadius;
      }
    else
      {
      radius[d] = ( (imageSize[d] - 1) / 2 );
      }
    }

  // Find the data-set boundary faces.
  FacesCalculator boundaryCalculator;
  FaceListType    faceList = boundaryCalculator(testImage, threadRegion, radius);

  // Support progress methods/callbacks.
  ProgressReporter progress( this, threadId, threadRegion.GetNumberOfPixels() );
  // Process the internal face and each of the boundary faces.
  for( FaceListIterator face = faceList.begin(); face != faceList.end(); ++face )
    {
    SmartIterator  test(radius, testImage, *face); // Iterate over test image.
    InputIterator  valid(validImage, *face);       // Iterate over valid image.
    OutputIterator out(outputPtr, *face);          // Iterate over output image.
    if( !test.GetNeedToUseBoundaryCondition() || !m_IgnoreBoundaryPixels )
      {
      test.OverrideBoundaryCondition(&nbc);
      for( valid.GoToBegin(), test.GoToBegin(), out.GoToBegin();
           !valid.IsAtEnd();
           ++valid, ++test, ++out )
        {
        // Get the current valid pixel.
        InputPixelType t = valid.Get();

        //  Assume a good match - so test center pixel first, for speed
        RealType difference = static_cast<RealType>( t ) - test.GetCenterPixel();
        if( NumericTraits<RealType>::IsNegative(difference) )
          {
          difference = -difference;
          }
        OutputPixelType minimumDifference = static_cast<OutputPixelType>( difference );

        // If center pixel isn't good enough, then test the neighborhood
        if( minimumDifference > m_DifferenceThreshold )
          {
          unsigned int neighborhoodSize = test.Size();
          // Find the closest-valued pixel in the neighborhood of the test
          // image.
          for( unsigned int i = 0; i < neighborhoodSize; ++i )
            {
            // Use the RealType for the difference to make sure we get the
            // sign.
            RealType differenceReal = static_cast<RealType>( t ) - test.GetPixel(i);
            if( NumericTraits<RealType>::IsNegative(differenceReal) )
              {
              differenceReal = -differenceReal;
              }
            OutputPixelType d = static_cast<OutputPixelType>( differenceReal );
            if( d < minimumDifference )
              {
              minimumDifference = d;
              if( minimumDifference <= m_DifferenceThreshold )
                {
                break;
                }
              }
            }
          }

        // Check if difference is above threshold.
        if( minimumDifference > m_DifferenceThreshold )
          {
          // Store the minimum difference value in the output image.
          out.Set(minimumDifference);

          // Update difference image statistics.
          m_ThreadDifferenceSum[threadId] += minimumDifference;
          m_ThreadNumberOfPixels[threadId]++;
          }
        else
          {
          // Difference is below threshold.
          out.Set(NumericTraits<OutputPixelType>::Zero);
          }

        // Update progress.
        progress.CompletedPixel();
        }
      }
    else
      {
      for( out.GoToBegin(); !out.IsAtEnd(); ++out )
        {
        out.Set(NumericTraits<OutputPixelType>::Zero);
        progress.CompletedPixel();
        }
      }
    }
}

// ----------------------------------------------------------------------------
template <class TInputImage, class TOutputImage>
void
LOCAL_DifferenceImageFilter<TInputImage, TOutputImage>
::AfterThreadedGenerateData()
{
  // Set statistics about difference image.
  ThreadIdType numberOfThreads = this->GetNumberOfThreads();

  for( ThreadIdType i = 0; i < numberOfThreads; ++i )
    {
    m_TotalDifference += m_ThreadDifferenceSum[i];
    m_NumberOfPixelsWithDifferences += m_ThreadNumberOfPixels[i];
    }

  // Get the total number of pixels processed in the region.
  // This is different from the m_TotalNumberOfPixels which
  // is the number of pixels that actually have differences
  // above the intensity threshold.
  OutputImageRegionType region = this->GetOutput()->GetRequestedRegion();
  AccumulateType        numberOfPixels = region.GetNumberOfPixels();

  // Calculate the mean difference.
  m_MeanDifference = m_TotalDifference / numberOfPixels;
}
} // end namespace itk

#endif
